Introduction

Easily recognizable traces of dramatic climatic variations make high mountain areas unique geotopes and “storytellers” about past as well as potential future climate change effects on landscape dynamics and living conditions in regions of rugged topography. Long-term observations of glaciers have provided convincing evidence of rapid global climatic change; the worldwide retreat of mountain glaciers during the 20th century was striking (Haeberli et al., 1999; Meier and Bahr, 1996). The apparent homogeneity of the signal at the secular time scale, however, contrasts with the large variability at local/regional scales and over time scales of years to decades (Letréguilly and Reynaud, 1990). Based on an increased theoretical and empirical knowledge of processes and feedback effects, it has been possible also to predict possible climate changes on a regional scale (IPCC, 2001; OcCC, 2002). Direct and dramatic ecological responses to this impending warming are expected (Peters and Lovejoy, 1992), in the form of feedbacks that could modify transfer rates of energy, water, and trace greenhouse gases at the earth's surface (Rosenberg et al., 1983). Landscapes may respond very noticeably and differentially to climate change as they integrate all ecological and historical factors (Theurillat et al., 1998). A key or “interface” function must thereby be attributed to soils. Soil-landscape patterns result from the integration of short- and long-term pedogeomorphic processes (Friedrich, 1996; Klingl, 1996). Despite numerous studies on the effects of climate warming on single processes, little is known about the reaction of a whole soil ecosystem (Rustad et al., 2001).

Many of the soil properties change continuously with time. The soil, however, can only be measured at a finite number of places and times, and any statement concerning the soil at other places or times involves prediction. Variation in soil is so complex that no description of it can be complete, and so prediction is inevitably uncertain (Heuvelink and Webster, 2001). Only appropriate observation will provide the basic knowledge necessary for assessing the development in reality, and statistically calibrated numerical spatial/gray-box models (rather than sophisticated deterministic models as applied in detailed scientific process studies) must help in the prediction of consequences and possible future scenarios. The results of Mendonça Santos et al. (2000) and Herrmann et al. (2001) demonstrate the utility of GIS technology for the facilitation of data-set management, for spatialization, analysis, and visualization.

Anthropogenic emission of greenhouse gases is predicted to increase the earth's average surface temperature during the next 50–100 years. For the area of interest (proglacial area of Morteratsch, Switzerland), a temperature increase of +1.6°C by the year 2050 and +3°C by the year 2100 can be assumed (according to the scenario A1B of IPCC, 2000). Therefore, additional areas will become ice-free and subject to weathering and soil formation. The most evident soil changes in the Alps will occur in proglacial areas where already-existing young soils will continuously develop and new soils will form due to the glacier retreat. The rate of the reactions that are of fundamental interest in the understanding of the soil system and its interaction with the environmental surrounding conditions has been deduced by a chronosequence in this proglacial area (Egli et al., 2006 [this issue]).

The main aim of this investigation was to predict future soil development for the next 100 years in the existing and new proglacial area associated with the retreat of the glacier Morteratsch and using the previous findings from chronosequences and statistical analyses (Egli et al., 2006 [this issue]).

Investigation Area

The studied soils lie within the proglacial area (1900–2150 m a.s.l., with a present-day area of 1.8 km²) of Morteratsch in the Upper Engadine (Switzerland). The investigation area is described in more detail in Egli et al. (2006 [this issue]).

Methods

Modeling of Soil Properties

Spatial modeling was performed with ArcGIS 8.3 (ESRI) with modules programmed in Visual Basic for Applications (VBA). Input data sets were the digital soil map, the glacial states, the digital terrain model (DTM; raster of 20 m) within the proglacial area, the digital elevation after the subtraction of the glacier surfaces (Biegger, 2004; see below), and simulated extents of the Morteratsch glacier for the years 2050 and 2100 (Biegger, 2004). The calculations were done raster-based (GRID; 20 m resolution).

Soil types were evaluated using the statistical trends of the individual soils as a function of time, landscape form, exposure, and slope according to Egli et al. (2006 [this issue]) and to Tables 1 and 2.

Table 1

Regressions between topographical features and time for the soil type Skeletic/Lithic Leptosol.

Table 2

Regressions between topographical features and time for the soil type Humi-Skeletic Leptosol.

The regressions are a mathematical expression of the probability of a certain soil type as a function of time and a given topographic feature. Equations 1–3 show how the probabilities of a certain soil type on a specific grid with certain topographic characteristics were calculated. The probabilities of a soil type were calculated as a function of the slope W_S, the exposure W_E, and shape of the landscape W_L. These probabilities are, furthermore, related to the time of the soil formation. The semi-empirical factors a, b, and c are derived from the regression equations.

1

2

3

The investigated factors are independent. To calculate the probability of a soil type (W_T) with a certain age at a specific grid, the probabilities according to Equations 1–3 have to be multiplied.

4

The minimum probability for the occurrence of a soil type has been determined in an iterative process (learning process; cf. also Behrens et al., 2005), i.e. an optimization of the modeled soil types with the soil map. The output of the model finally produces values ranging from 0 to 1. The closer the values are to 1, the more likely is the occurrence of the considered soil type. The classification returning the highest prediction accuracy is used for the final prediction.

The soil properties are strongly dependent on the soil type, shape of the landscape, topography, and age. A more detailed description of the modeled soil properties is given in Table 3. The modeling of the soil properties, consequently, was carried out heuristically. The modeled soil types served as a base that was related to the soil age and topography to derive the individual soil properties. The calculation was according to the entity-relationship principle (Klingl, 1996). The implementation was done by a query of the databank using Boolean functions such as “and,” “or,” “not” as well as logical operators such as “lower than,” “greater than,” “equal,” “not equal” to get the corresponding properties (Klingl, 1996). The pH values in the topsoil (surface soil layer) and subsoil, the soil depth relevant for plant growth, and the skeleton content in the top- and subsoil could be modeled with this procedure (see Table 4, which lists the model structure for the calculation of the specific soil properties).

Table 3

Modeled soil properties.

Table 4

Modeling of soil properties (entity relationships containing Boolean and logical functions) derived from soil types, soil age, and topography.

Soil cartography and classification (also of soil properties) was made according to the FAL system (Brunner et al., 1997). Soil types are translated into the WRB (FAO, 1998) system. The dataflow is given in Figure 1. Some specific properties were modeled according to the FAO-UNESCO system (1990) as the WRB system does not fully meet the requirements for a soil classification in Alpine areas; e.g. the soil type “Ranker” indicates an intermediate stage between Humi-Skeletic Leptosols and Dystric Cambisols (endoskeletic).

Figure 1

Dataflow used for the modeling of soil types (a) and properties (b).

¹Soil depth: relevant for plant growth.

Modeling the Surface Area Underneath the Glacier

The present topographical situation will change with a further retreat of the Morteratsch glacier. To model soil properties in a future proglacial area (that is presently covered by ice), it is essential to have or infer information about the surface underneath the glacier because one important model input is the topography (see above). To calculate this topography the glacier was regarded as the sum of several surfaces, each of them representing an individual branch. A cubic Non-Uniform Rational B-Spline (NURBS) (Piegl and Tiller, 1997) curve was applied to parameterize the course of a branch in flow direction (Biegger, 2004). A NURBS surface S(u,v) of degree p in the u direction and of degree q in the v direction is a piecewise rational function

5

where n+1 (m+1) are the number of control points in u (v) direction, {P_i,j} form a bidirectional control net arranged in a rectangular fashion, and {R_i,p} ({R_j,q}) are the rational functions in the u (v) direction. Typically, a mountain glacier can be thought of as a combination of smaller parts called branches where the number of branches n_b varies for each glacier. The upper boundary of the longitudinal profile is defined by the parametric flow line C_f. The flow line runs from the source point P_s to the terminal point P_t. The vertical dimensions of the longitudinal profile are calculated using steady-state conditions as proposed by Haeberli et al. (1999). The glacier thickness along the flow-line C_f can be expressed as a combination of the ice cap factor k_x and the relief factor r_x. The maximum thickness d_max of the profile would result at the equilibrium line (ELA = equilibrium line altitude of glaciers). The ice cap factor can be written as (Biegger, 2004)

6

7

where x is an index indicating the position on the flow-line relative to P_s, l_x represents the length from P_s to the surface point P_x at x, l_ELA is the distance from the point P_s to P_ELA, and L = $\overline {C_f } $ is the glacier length measured along the sloping surface using an empirical relation between glacier thickness and glacier width. A continuous description of the branch surface is achieved by approximating the parabolic cross sections with a NURBS surface S(u,v) of degree 2 in the transverse direction and of degree 3 in the flow direction (Biegger, 2004). To geometrically model the branch surface, the longitudinal profile had to be extended with n_x cross sections. Although a glacier may consist of more than one branch, each branch had to be treated separately. Each branch runs from the source point P_s to the terminal point P_t. In order to compute the glacier bed, the surface S(u,v) was discretized by projecting each grid point of the original DEM onto the branch surface and subtracting the thickness from the DEM at each grid point. For the deformation of a glacier due to climate change, each branch was considered individually and its reaction was parameterized as proposed by Haeberli and Hoelze (1995). Climate change was introduced numerically by changing the mass balance of the glacier. Given a step-wise temperature change ΔT_Air assuming a certain climate scenario, the change Δb of the mass balance could be calculated by

8

where ΔELA/ΔT describes the vertical increase of ELA per 1°C warming and integrates the change of all climate parameters (humidity, radiance, accumulation and air temperature, feedback effects; Kuhn, 1990). Experimental studies have shown a good accordance with historical measurements applying a value of 160 m/°C for the Morteratsch glacier. The experimentally evaluated ratio of Δb/Δh was 0.68 (Biegger, 2004).

Results

The spatial distribution of the simulated glacier thickness is visualized in Figure 2. It can be seen that there are two regions of increased ice thickness. Both regions are located in a plane just below steep slopes where the thickness is less than 100 m. The applied spatial sampling was 20 m according to the horizontal resolution of the DTM. The surface area with and without the Morteratsch glacier (glacier bed) is plotted in Figure 3. The future evolution of the Morteratsch glacier distinctly depends on the climate scenario. The sensitivity study by Gyalistras (2000) for estimating the future trends of air temperature and precipitation in the Alps showed that air temperature will possibly rise more than the global average. A probable climate scenario could be the one according to the scenario A1B of IPCC (2000) with a temperature increase of +1.6°C by the year 2050 and +3°C by the year 2100. The impact of the applied climate scenario on the Morteratsch glacier is visualized in Figure 4. The most obvious changes can be recognized at the tongue, whereas the accumulation area remains more or less unaffected. In the year 2100 the branches formerly contributing to the Morteratsch glacier will become isolated. Within the next 100 years a much greater area than the present-day proglacial area will become ice-free where soils will be able to form. According to the surface area calculations, new proglacial lakes are not supposed to be formed in that time span.

Figure 2

Contour plot of the thickness of the Morteratsch glacier (1985; Biegger, 2004).

Figure 3

Present-day surfaces with (a) and without (b) the Morteratsch glacier (glacier bed). DEM25 reproduced by permission of swisstopo (BA067583).

Figure 4

Future development of the Morteratsch glacier assuming a temperature increase of +1.6°C until the year 2050 and +3°C until 2100.

DEM25 reproduced by permission of swisstopo (BA067583).

The basis for soil modeling included the time since deglaciation and topographic features. The slope, shape, and exposure of the landscape were calculated using the DTM25 (with 20 m resolution) and the modeled surfaces according to Biegger (2004). Additionally, the water net and the flow direction (modeled surfaces) were taken into account. The soil model calculates in the present proglacial area larger areas having Skeletic/Lithic Leptosols after about 30 to 50 years. Already after 90–100 years of soil evolution a transition (during a time span of 20–40 years) of Skeletic/Lithic Leptosols into Humi-Skeletic Leptosols is modeled which agrees well with the soil map. Along the main river system, the distribution of the Fluvisols generally agrees well with reality. A lower agreement was, however, found along small river branches. A comparison between the modeled area with soil cover and mapped soil area gave an agreement of about 74% (see also Fig. 5). Problems arose especially where the sediment bed of the glacier was extremely thin or even absent.

Figure 5

Optical comparison between mapped (left) and modeled (right) soils in the proglacial area. The black lines indicate schematically retreat-isochrones (retreat of the glacier since 1850) until 1997. The glacier state is for the year 2000. DEM25 reproduced by permission of swisstopo (BA067583).

Several soil characteristics have been modeled and compared with the soil map. Among the modeled characteristics are the pH value (CaCl₂) in the topsoil (uppermost, surface layer) and in the subsoil, skeleton content in the top- and subsoil, and the soil depth relevant for plant growth (according to the FAL system: soil depth = soil volume – skeleton volume – groundwater volume; result is related to depth instead of volume). In the correctly modeled soil area, the agreement of the modeled characteristics with the mapped ones was in all cases very high (Table 5); only very small differences to the soil mapping could be discerned.

Table 5

Comparison between modeled and mapped soil properties (values refer to the correct modeled soil area).

The stream network of the future proglacial area had to be calculated using the modeled surface area under the glacier. Together with the existing stream network the future situation was derived. The stream network was a basis for the calculation of the distribution of Fluvisols.

During the first 25 years of soil formation only patches having a significant soil cover could be found or are modeled. Thus, a certain area in front of the glacier tongue remains without any soil cover. The area with soil cover then steadily increases with Skeletic/Lithic Leptosols that appear first and transform later into Humi-Skeletic Leptosols (Fig. 6). The older Humi-Skeletic Leptosols show signs of an initial B horizon formation (here called Ranker; see Fig. 6). This soil type will increasingly appear in the year 2050. Ranker will first be formed in rather flat areas with a soil age of >150 years. The eastern part of the investigation sites should develop slightly more quickly due to having more north-facing slopes. After 200 years (maximum age in the year 2050) of soil development, nearly 96% of the glacier bed will be covered by a (thin) soil layer. A major problem in soil modeling is the consideration of the lateral moraines. They are usually very unstable and disturb soil formation. The periphery of the lateral moraines has to be considered as an area more for potential than for effective soil formation. Ranker will be the dominant soil type in the front of the proglacial area at the end of the 21st century (Fig. 6). In the upper part of the proglacial area mostly Skeletic/Lithic Leptosols and Humi-Skeletic Leptosols will be found. In flat parts close to the main river, additional Fluvisols should develop. A considerable part of the upper proglacial area does not have any soil cover. Lithic/Skeletic to Humi-Skeletic Leptosols are modeled on the young lateral moraines. Any soil development there will depend on the stability of these moraines ( = potential areas of soil formation).

Figure 6

Soil development in the proglacial area in the next 150 years. The predictions are based on the assumption that temperature will increase +1.6°C until the year 2050 and +3°C until 2100. Modeling for the year 2150 is based on the assumption of a continuous retreat (as no prediction of the equilibrium for that time was available). DEM25 reproduced by permission of swisstopo (BA067583).

There was no model for the glacier tongue states available for the year 2150. Assuming that the glaciers will not have readvanced to the original position of the year 1850 ( = lowest part of the proglacial area), which is highly improbable, soil evolution can be calculated at least for the lower part of the proglacial area. If we assume that the glaciers will retreat also in the 22nd century, then a soil development can also be modeled for the upper part of the proglacial area (Fig. 6). After 300 years of soil evolution (lower part of the proglacial area, year 2150), Dystric Cambisols (endoskeletic) will probably appear. Because this soil type has not been found in the present-day proglacial area (after 150 years of soil evolution), its first appearance was derived from chronosequences in the central Alps (Egli et al., 2001). This makes the prediction for the year 2150 even more speculative and was, therefore, only done for the soil types.

Discussion and Conclusion

Soil evolution will relatively rapidly follow the glacier retreat. However, it takes about 25 years until significant signs of soil formation will be observable. Within the next 100–150 years the soils in the present-day proglacial area will thicken; the skeleton content will be reduced by being subjected to physical and chemical weathering (and also because of the accumulation of organic matter; Fig. 7), will continuously acidify (Fig. 8), and will progressively melanize (darken) with time and distance from the ice. Soil types will grade finally into Rankers (Humi-Skeletic Leptosols with traces of a B horizon) and most probably not before the year 2100 into Dystric Cambisols (endoskeletic). Soil evolution, although quick in the Alps, seems to be slower when compared to similar sites in Norway where Regosols were observed already after 47 years of deglaciation which then graded into Brunisols (Cambisols) in less than 120 years (Mellor, 1985, 1986). Frost disturbance as well as slope instabilities can, however, impede or delay such a progression of soil types (Haugland, 2004). This can lead to highly patterned features of the spatial soil distribution associated with an abrupt threshold between genesis and stabilization. As a consequence, modeling of soils on a small scale in existing and future proglacial areas can hardly be very precise.

Figure 7

Modeled and predicted skeleton content (volume-%) in the topsoil for the next 100 years. DEM25 reproduced by permission of swisstopo (BA067583).

Figure 8

Modeled and predicted pH values in the topsoil for the next 100 years; strongly acid corresponds to pH (CaCl₂) 3.3–4.2, acid to 4.3–5.0, and weakly acid to 5.1–6.1. DEM25 reproduced by permission of swisstopo (BA067583).

Modeling alpine soil properties can be done using several methods. Huber (1994) made a process-oriented modeling of soils with GIS in the Bavarian Alps using topographic indications, comparisons with existing soil maps, and the formulation of geo-ecological interrelationships based on the principles of Leser and Klink (1988). A further possibility is the derivation of soil relevant data using remote-sensing techniques (e.g. Mikhaylov, 1990; Gauthier and Tabbagh, 1994). If enough point data are available, then geostatistical principles can also be used for predicting soil properties (e.g. Ahn et al., 1999; Grunwald et al., 2000; Zhu et al., 2001; Lark, 2003). The applicability to larger regions of detailed soil surveys in small reference areas was successfully assessed by Lagacherie and Voltz (2000). The mapping method consisted of a soil classification in a reference area. The probabilities of occurrence of soil classes at a site were used to predict soil properties at unvisited sites using digital data sets such as the DTM, geology, and hydrology. Gessler et al. (2000) were able to account for between 52 and 88% of soil property variance using easily computed terrain variables such as slope and flow accumulation. Empirical models developed from DTM were also successful in predicting horizon depths of the topsoil (Thompson et al., 2001). Promising results in mapping soil properties are also obtained from fuzzy functions (Hannemann, 2005) or neuronal networks (Behrens et al., 2005).

One drawback of all these mapping techniques is the absence of the time scale that was important for the case study presented. Chronosequences and therefore the relationship of soil properties to time as a function of topographic properties are vital in making any (4D) predictions of soil evolution in proglacial areas. The statistically and probabilistically based model, however, also has its weaknesses. The modeling of the present-day situation of soil distribution gave “only” an agreement of about 73% with the soil map, although very significant correlations of the individual soil types with topographic features and time could be found. One major problem occurs because the sediment bed of the glacier was partially extremely thin or even absent. In such a case, soil modeling failed. Another problem is related to the lateral moraines that may be unstable and consequently hinder soil development. Erosion or small debris flows from such moraines occur randomly and were not really predictable. Finally, the interaction of geomorphodynamics and soil formation makes soil modeling difficult in high Alpine environments as described in Egli et al. (2006 [this issue]).

Nevertheless, the modeling of soils will, for instance, allow a more precise prediction of vegetation changes and thus enable a more comprehensive prediction of future landscape evolution in a rapidly changing high Alpine environment.